arXiv:1509.00322v2 [cond-mat.mes-hall] 3 Oct 2015 


uiaiL 


Multidomain switching in the ferroelectric nanodots 

PiERRE-WiLLiAM Martelli^, Seraphin M. Mefire^ and Igor A. Luk’yanchuk^’^ 

^ University of Lorraine, lECL, CNRS UMR 7502, 54506 Vandoeuvre-les-Naney Cedex, Franee 
^ University of Pieardie, Laboratory of Condensed Matter Physies, 80039 Amiens Cedex 1, Franee 
^ L. D. Landau Institute for Theoretieal Physies, Moseow, Russia 


PACS 77.80. Dj - Domain structure; hysteresis 
PACS 64.60. an - Finite-size systems 

PACS 02.60. Cb - Numerical simulation; solution of equations 

Abstract -Controlling the polarization switching in the ferroelectric nanocrystals, nanowires 
and nanodots has an inherent specificity related to the emergence of depolarization field that is 
associated with the spontaneous polarization. This field splits the finite-size ferroelectric sample 
into polarization domains. Here, based on 3D numerical simulations, we study the formation of 
180° polarization domains in a nanoplatelet, made of uniaxial ferroelectric material, and show that 
in addition to the polarized monodomain state, the multidomain structures, notably of stripe and 
cylindrical shapes, can arise and compete during the switching process. The multibit switching 
protocol between these configurations may be realized by temperature and field variations. 


A fundamental property of ferroelectrics is the interplay 
between the spontaneous polarization, P, appearing due 
to the symmetry-breaking off-center displacement of po¬ 
lar ions, and the long-range depolarization electric field 
caused by the same polarization. The depolarization field 
is induced by depolarizing charges distributed with vol¬ 
ume density p = divP in regions where the polarization is 
nonuniform and with surface charge density a = Pn in the 
near-surface layer of polarization termination points (here, 
n = {ux^ny^Uz) is the unit vector, normal to the sample 
surface and directed outside the sample). As illustrated in 
Fig.[TJi for the finite-size sample, the depolarization field 
is distributed in the inner space and in the surrounding 
outer space that costs an additional electrostatic energy 
and impedes the formation of the ferroelectric state. 

The tendency to reduce the unfavorable action of the 
depolarization field was noted by Landau and Kittel more 
than 60 years ago (TJ|^ as able to lead to a polarization 
domain patterning of the sample. As shown in Figs. Etc, 
the formation of the oppositely polarized 180° domains al¬ 
ternates the associated depolarizing charge at the surface. 
This confines the depolarization field closer to the surface, 
diminishing its energy. The price for the gain however is 
the additional cost of domain wall (DW) creation. Hence, 
the domain form and size are the result of a tiny balance 
between geometry, electrostatic and ferroelectric proper¬ 
ties of the system. Figs, [^c exemplify the competition 
between the stripe and cylindric domain patterns. The 


situation becomes even more delicate if an external field 
is applied. The interaction with the field changes the do¬ 
main configuration, favoring the ” up”-oriented domains 
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Fig. 1: Distribution of the depolarization electric field in the 
ferroelectric nadodot (ND) with the polarized monodomain 
state, MDS (a), the two-band state, TBS (b), and the 
concentric-cylindrical state, CCS (c). Panel (d) shows the 
multibit switching between these states, denoted by quantum 
numbers 0, ±1 and ±2. The red and blue colors correspond to 
the ”up” and ’’down” polarization orientations. 
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with polarization parallel to the field. The growth of the 
up-polarized domains continues up to complete poling of 
the sample to the monodomain state. 

The formation of Landau-Kittel domains having the 
structure of regular periodic stripes in thin ferroelectric 
films was studied over the past decade both experimentally 
|4l|^ and theoretically [8-11 


The period of the domain 27 28 


2d ~ 2^3.53 (e/£c) 


films 18 confirm this hypothesis. 


19 


21 


high-density memory-storage units 13 . The experimental 
and theoretical ^2 2b efforts were mostly concen- 


of the near-surface dead layer 26 


In this letter we study how the Landau-Kittel struc¬ 
ture of 180° polarization domains is formed in a finite-size 
sample. We show that the confining electrostatic effects 
result in the various domain structures. Field and temper¬ 
ature applications permit to realize the controllable multi¬ 
bit switching between them, hence increasing the volume 
of the writable information per ND. 

To disjoint the depolarization effects, that are primary 
for the domain formation, from the secondary lattice- 
deformation effects, we consider the uniaxial ferroelectric 
material for which the ferroelastic coupling is small. To be 


specific, we select the Sodium Nitrite, NaN02, as model 
material. It has the orthorhombic symmetry, hence only 
two, ”up” and ’’down”, orientations of the spontaneous po¬ 
larization Pq 12 fiC/crn^ with respect to the c-axis. The 
anisotropic tensor of dielectric constants has the principal 
values Sa — '^0, Cb — b and 5c — 12 at low temperatures 
providing the preferred orientation of DW in the 


pattern was found to scale according to the Landau-Kittel 
square root law as js 12 


6c-plane. 

The ferroelectric transition in the bulk NaN02 crystal 
occurs at Tc ^ 433 K, the Curie constant, C, being equal 
The characteristic value of the depolariza- 


(1) to 5000 K 27 


where h is the film thickness, Co — ^ atomic-scale 

coherence length, 5 = 5p -h ^.nd Sa are the in¬ 

trinsic longitudinal and transversal (to polarization and to 
stripes) dielectric constants of the ferroelectric state and 
Sp is the dielectric constant of the surrounding paraelec- 
tric environment. When the electric field, Ea^ is applied, 
the theoretical calculations predict the stretching of up- 
oriented domains and the contraction of down-oriented do¬ 
mains. Interestingly, this induces the oppositely-oriented 
coarse-grained depolarization field that finally leads to the 
negative permittivity of the ferroelectric layer pp5] . The 
behavior of the domain structure at higher fields is lit¬ 
tle studied although one can guess that the oppositely- 
oriented domain stripes will transform to the vanishing 
domain droplets before the complete polarization of the 
system, similarly as it happens in ferromagnetic samples 
16, 17 . Recent numerical calculations for ferroelectric 


Much less is known concerning the properties of the 
Landau-Kittel 180° domains in finite-scale nanodot (ND) 
samples, with sizes comparable to the predicted domain 
width d. Meanwhile, it is the application of the NDs 
that is considered to be very promising in the emerging 
ferroelectric-based nanoelectronics for the realization of 


trated on the pseudo-cubic perovskite ferroelectrics with 
eight possible polarization orientations. The diversity of 
domains and of polarization-vortex patterns was discov¬ 
ered, but the systematic study was partially impeded by 
the complexity of the system. The electrostatic depolar¬ 
ization effects are strongly coupled with ferroelastic ones, 
produced by the polarization rotations and by the strain 


tion field is estimated as Eq = Pq/sq 13.5 x lO'^kV/cm 
(here Sq is the vacuum permittivity). We assume also that 
the crystal is cut in the form of cylindrical nanoplatelet 
of radius r 8nm and of thickness h 1.7 nm, the 
spontaneous polarization being directed along the cylinder 
axis. The ND is embedded into the dielectric matrix with 
Sp 90. The high value of Cp is required to reduce the de¬ 
polarization energy. Otherwise the depolarization effects 
will kill the ferroelectric phase, even in the multidomain 
state. With such a selection of parameters, the charac¬ 
teristic domain width, calculated according to Eq. 0, is 
d 7.1 nm that is commensurate with the ND radius, r. 

The principal results are illustrated in Fig. At zero 
applied field the ND can stay either in the polarized mon¬ 
odomain state (MDS), Fig. [^, or in one of the mul¬ 
tidomain states. The competition occurs mostly between 
the two-band state (TBS), Fig. [^, and the concentric- 
cylindrical state (CCS), Fig. [^, although other config¬ 
urations are also possible. The average polarization of 
the TBS, P, is zero whereas the CCS and the MDS have 
the nonzero polarizations, either “up”-oriented or ”down”- 
oriented. We ascribe the appropriate quantum numbers, 
= 0, ±1 and ±2, to the TBS, CCS and MDS respec¬ 
tively, where the sign reflects the orientation of P. 

The field application allows to jump among the TBS, 
CCS and MDS as sketched in Fig. [^. Usually, the TBS 
arises at zero-field cooling. When the field, Ea > 0, is 
applied the TBS stays reversibly-stable until it jumps to 
the up-oriented MDS. Further field increase and reversal 
cycling produce the hysteresis curve that, depending on 
the field variation protocol, realizes the four-bit switching 
between the CCS with N = ±1 and the MDS with N = 
±2. Although the initial TBS with = 0 is no longer 
accessible, it can again be achieved by ^dhermal reboot^’ 
consisting in the heating of the system to the paraelectric 
state with subsequent zero-field cooling. The interplay 
between the TBS and the CCS during the ND poling is 
a legacy of the already mentioned field-induced transition 
between stripe and droplet domain configurations in an 
infinite system. 

An analogous multiple-hysteresis phenomenon was re¬ 
cently observed in mesoparticles of superconducting lead 
and explained by an irreversible decay of magnetic vor¬ 
tex droplets 29 . The effect was shown to be similar to 


the Rayleigh fragmentation of a charged liquid due to the 
competition between the long-range Coulomb charge re- 
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pulsion and the short-range molecular attraction 30 . The 


buckling, faceting, or even disintegration of ferroelectric 
domains also can be driven by the long-range repulsion 
between DWs, associated with the electric fringing fields 
of depolarization charges 31 


To give the further insight on the switching process we 
describe the methodology of calculations. We consider 
the nonuniform Ginzburg-Landau equation for the order 
parameter, spontaneous polarization P, coupled with the 
electrostatic equation for the electric potential ipf, in flj, 
the interior of the ND, 


10 


11 


i - Clbdy - P = -soKcdziff, 

So {sad‘^ -h Shdy + Sicd‘1^ iff = dzP. (2) 


Here, t = (T — Tc)/Tc is the reduced temperature, /^c = 
C /Tc — 11.5 is the displacive parameter |^, Sic ^ 1 is 
the small nonpolar contribution to the Curie-Weiss diver¬ 
gency of 5c, such as Scit) = Sic + /^c/^ at T > Tc. The 
coordinates x, y and 2 : are selected along the axes a, b and 
c, the corresponding coherence lengths C06 ^Oc 
roughly equal to 1 nm. To catch the generic features of the 
switching we use the simplistic cubic form of the nonlin¬ 
ear term, digressing from the more specific details of the 
transition. The parameter Pq is selected as spontaneous 
polarization at low temperatures. The electric field and 
the electric displacement are given by Ej = and 

Y>f = 5oEj + P/, where the components of the polariza¬ 
tion vector, Pj, are written as: Pf^x = ^0 (^a ~ Pf,x^ 

Pf,y — '^0 ('^6 1 ) Pf,y and Pf^z — '^0 ('^ic 1) Pf,z T P • 

The relations § are completed by the Poisson equa¬ 
tion = 0 for the electric potential (pp^ outside the 

ND. The corresponding electric field and displacement are 
expressed as Ep = —Vifp and Dp = 5o5pEp. 

Electrostatic rules imply the continuity of the poten¬ 
tial, iff = ipp, and of the normal component of the 
electric displacement, Djn = Dpii. In our notation, 
the latter constraint is written as SoSicdzP^f — P = 
^o^pdz^p for the upper and lower ND surfaces and as 

{saTixdx + SbUydy) ipf = Sp {rixdx + riydy) (Pp for the lat¬ 
eral surface. For the polarization we took the free bound- 
ary condition {il^nxdx + ilb^ydy + il^nzdz) -P = 0 at the 
sample surface. 

The calculation setup is shown in Fig. To model 
the outer potential distribution, the ND was embedded 
into a bigger space, ^2, with sufficiently large volume to 
avoid to disturb the emergent fringing field. In practice the 
cylinder with radius R cs 12 nm and thickness H cs 17 nm 
was taken, the independence of the obtained results on the 
cylinder sizes was checked at each calculation stage. The 
potential difference A(pp = U between the top and bottom 
electrodes was applied to induce the field Ea = U/H^ the 
variation of cpp at the lateral surface being assumed to be 
linear, 7 ^p{z) = {z/H)U. 

We deal then with the boundary value problem satisfied 
by P and (p^ where the electric potential in ft is denoted 


by (^, as being or Pp, namely inside or outside the 
ND. The approach for solving numerically this 3D prob¬ 
lem is based on a finite element method, implemented as 
Fortran 90 homemade code (see 33 for details). First, 


we introduced a variational formulation of the problem, 
for which the variational unknowns, representing P and 
p, are found in the Sobolev spaces H^{ftf) and H^{ft) 
respectively 34 . Then, we discretized this formulation 


by making use of a mesh of ft (the closure of ft) and 
of the Lagrange finite elements of the first order. This 
mesh consists of a collection of tetrahedra, obtained from 
a usual process of triangulation, and is fine enough near 
the boundaries of ftf and ft. The discrete regions associ¬ 
ated with Qf and Q are polyhedral. In a correlative way, 
the discrete region associated with the boundary of ftf 
is entirely made up of faces of tetrahedra. Each of these 
faces is common both to a tetrahedron of the discrete re¬ 
gion associated with ftf and another one of the discrete 
region associated with ftp = ft\ftf. By dealing with a 
mesh size equal to 0.47nm for ftf and to 0.61 nm for ftp, 
we were led to a square nonlinear system of 485621 scalar 
equations. Let us mention that this system is not subject 
to a uniqueness of solution, as it is the case for the as¬ 
sociated discrete and continuous variational formulations. 
Finally, we solved the system with the help of an inexact 
Newton method 35 , combined at each iteration with the 


Generalized Minimal RESidual (GMRES) algorithm 36 


a preconditioner, based on an incomplete LU factorization 
(see, e.g., [^), was incorporated into this algorithm. 

We systematically dealt with the initialization of the 
inexact Newton method either with a randomly generated 
datum or with a solution profile obtained at the previous 
calculations. The behavior of the system as a function of 
the temperature or of the field was followed with the step 
of at most 10“^ of the running interval whereas the finer 
step of 10“^ was used in the points of the brutal solution 
changes. 

(§ were also studied in 


Eqs. 


38 for the ellipsoidal 


ND of triglicine sulfate (TGS), having similar dielectric 
properties. In particular the temperature-induced GGS- 
MDS transition was observed. However, only the axially- 
symmetric domain profiles (that are not always the most 



Fig. 2: Geometrical parameters of the system. 
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Fig. 3: (a)-(f) Different (meta)-stable polarization profiles, ob¬ 
tained from random initial polarization distributions. Lower 
panel shows the polarization distribution in the CCS for the 
equatorial (g) and vertical (h) sections of the ND. 
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Temperature, T/T, 


Fig. 4: Temperature evolution of the average polarization of 
the ND for different domain configurations. Dots denote the 
most stable state. The polarization is given in units of the 
zero-temperature bulk spontaneous polarization, Pq, the tem¬ 
perature is given in units of the bulk critical temperature, Tc. 



stable) were obtained. For instance, the steady TBS was 
overlooked. This does not permit to see all the diversity 
of the possible domain structures and to investigate the 
field-induced switching between them, that is the subject 
of the current article. 

In what follows, we start from study of temperature evo¬ 
lution of the system and show that the different metastahle 
domain states can be formed at the same temperature; 
some of them are presented in Fig.|^ To catch all the pos¬ 
sible states, we simulated the zero-field thermal quenching 
to low temperatures. Random polarization and random 
potential were taken at the initialization stage. Then, the 
final converging state was registered. Repeating the nu¬ 
merical experiment with different initial distributions, we 
found that, in addition to the already mentioned MDS, 
TBS and CCS (Figs. Ht-c, respectively), other profiles, like 
the asymmetric-droplet state (Fig. §1) and the three-band 
state (Fig. [^), can arise at T = 0. Several other domain 
patterns were also seen, but on a less regular basis. In the 
TBS, the polarizations of up- and down-oriented domains 
are compensated, giving P = 0. In other states the av¬ 
erage polarization is nonzero; the latter being maximal in 
the MDS. 

To follow the temperature evolution of the observed do¬ 
main states we performed the adiabatic heating when the 
next-on-heating state was obtained from the previous one, 
taking the latter as the initial configuration. The average 
spontaneous polarization of each of the discovered states 
decreases on heating (see Fig. [^, each state having its 
own temperature range of existence above which it irre¬ 
versibly jumps to another one, more enduring state. The 
corresponding transition temperature was identified as the 
temperature of superheating. Remarkably, no jump cor¬ 
responding to a supereooling instability was observed on 
adiabatic cooling. 

The MDS stays stable on heating till the temperature 
O.lSTc, above which it suddenly jumps to the CCS. The 
asymmetric-droplet state also transforms to the CCS by 


almost continuous coalescence of the edge-connecting neck 
but at higher temperature, 0.62Tc. The average polariza¬ 
tion of the CCS itself gradually decreases and at 0.63Tc 
the system transforms to the TBS. 

The TBS is the most enduring state. Its average po¬ 
larization is always zero, whereas the maximal domain 
polarization also decreases on heating and vanishes at 
To = 0.68Tc. This temperature can be considered as 
the ferroelectric transition temperature of the ND. As 
concerns the three-band state, it remains metastable till 
0.37Tc. Then, it transforms to the new diamond-like state 


(Fig-i) 


with P = 0, which at 0.65Tc drops to the TBS. 
We identified the most stable domain configuration 
at each temperature, calculating the free energy, F = 
— [ll|, for each state (see the dot marks in 

Fig-i- The TBS remains the most stable on cooling from 
To to 0.17Tc and then, the avail in energy passes to the 
MDS (see dashed line in Fig. [^. This thermodynamie 
transition temperature is slightly lower than the given 
above superheating temperature for the MDS, O.lSTc. 

Figs. illustrate the polarization distribution inside 
the CCS at low temperatures. The width of the DW is of 
the order of 4-6 ^o which is larger than the typical atom¬ 
istic DW width in bulk materials. The DW has Ising- 
like structure when the polarization changes in amplitude, 
across the wall, but stays almost parallel to the z-axis. At 
higher temperatures the DW profile becomes even softer 
10 11 . There is also no substantial change in the polar¬ 


ization value when passing across the ND from bottom to 
top. This makes a difference with the high-5 perovskite 
ferroelectrics, where the polarization flux-closure at ter¬ 
mination points of the DW has been observed [9 [mill- 
The determining difference is in the opposite ratio Sp/Sc 
that substantially reduces the depolarization effects in the 
embedded ND of NaN02. Another striking feature we 
observed is the deviation of DWs from the most prefer¬ 
able 6c-plane, clearly seen in the diamont-like structure 
(Fig- §), and also the DWs buckling in the three-band 
state and, especially, in the asymmetric-droplet state. We 
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Fig. 5: Hysteresis curves for different domain states at T = O.lTc. The average polarization of the ND is given in units of the 
spontaneous polarization, Po at T = 0, the applied electric field is given in units of the bulk depolarization field Eq = Pq/sq. 


believe that this effect is provided by the already men¬ 
tioned interplay between the positive DW tension energy 
and the negative nonlocal electrostatic energy of the fring¬ 


ing field 29 


The behavior of the multidomain structures in external 
field was also studied by a similar method of adiabatic field 
variation. Fig. (which is the extended version of Fig.[^) 
presents the hysteresis loops, P{Ea), for various domain 
states at T = O.lTc. Monitoring the domain evolution 
in the applied field gives more insight into dynamics of 
the switching process. Starting from the zero-field cooled 
TBS, having P = 0, and gradually increasing the field, 
we follow the central (blue) branch of the hysteresis loop. 
We notice that the DW in the TBS starts to bend, to 
diminish the unfavorable polarization. At some threshold 
field it irreversibly escapes from the sample to form the up- 
polarized MDS, which is the most stable state at > 0. 
The MDS remains with further field increase, forming the 
upper (red) branch of the hysteresis loop. Backward field 
decrease keeps the MDS till Ea = 0 and even further, for 
small negative Ea^ when the up-polarized MDS turns to 
be metastable with respect to the down-polarized MDS. 

Polarization distribution in the MDS is nonuniform. 
Feedback action of the depolarization field reduces the 
spontaneous polarization closer to the centers of the 
ND circular plane surfaces, the effect known as Landau- 
Lifshitz branching . Such a polarization central sagging 
becomes especially pronounced when the field downturns 
below zero. Finally, the polarization suddenly flips at the 
central part of the ND with formation of the (metastable) 
CCS. Notably, the average polarization of the just formed 
CCS is smaller than that in the original MDS but is still 
up-oriented, that is against the applied field. 

Further evolution of the system, shown by the brown 
line, depends on the protocol of the field cycling. While 


the down-oriented field continues to grow, the relative vol¬ 
ume of the internal negatively-polarized cylindrical do¬ 
main also increases. At some critical field, that is nothing 
but the coercive field of the system, the internal domain 
suddenly expands to the whole volume with the formation 
of the stable down-polarized MDS. This accomplishes the 
semi-loop of the global hysteresis, that can be completed 
by its centrally-symmetric counterpart under the rear¬ 
ward field sweep. If, however, right after the MDS^CCS 
switching, the running of the applied field reverses from 
the decrease to the increase, the internal domain dimin¬ 
ishes and, at some threshold > 0, suddenly vanishes to 
restore the original MDS. This closes the local hysteresis 
loop that can be observed between the MDS and CDS at 
small field oscillations. 

The evolution of the three-band and asymmetric-droplet 
domain states in the applied field (green and magenta 
lines in Fig.respectively) can be characterized as a ’’re¬ 
pulsive cycling”. Once these configurations are formed 
by the thermal quenching, they can live durably in their 
metastable state and be locally-stable under small field 
variation. However, in larger fields, they become unstable 
with respect to transition to the more regular TBS, CCS 
or MDS. As soon as this happens, there is no way to re¬ 
store the original state back neither via the temperature 
nor via the field variation. They disappear very soon af¬ 
ter the regular high-amplitude field cycling. For instance, 
one of the bent DWs of the three-band state escapes on 
the field variation and the system jumps to the TBS. Note 
the unusual crossing of the green and blue branches of the 
hysteresis curve, implying that P is not a unique function 
of Ea but depends also on the domain configuration. 

The behavior of the asymmetric-droplet state is even 
more diverse. Depending on the direction of the field 
variation, it can either join the CCS, similar as in the 
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temperature diagram in Fig. or jump to the TBS or 
to the MDS. The first case occurs when approaching the 
coercive field in which the system drops immediately to 
the oppositely-polarized MDS state. In the second case 
the accidental bifurcation degeneracy with respect to se¬ 
lection of the switching trajectory was observed. Within 
the temperature-induced fluctuation incertitude both the 
transitions to the TBS or to the MDS can occur (see 
dashed magenta line in Fig [^. The interesting feature 
of the transition to the TBS is the lowering of P when the 
applied field increases. 

The observed variety of nanoscale domain states with 
their irreversible history-dependent evolution can eluci¬ 
date the dynamics of the relaxor materials that, at certain 
conditions, can be viewed as a set of interacting polariza¬ 
tion nanoclusters with intrinsic domain structure, similar 
to the ferroelectric NDs. Parameters of our model system 
were selected in a special way to highlight the controlled 
five-bit switching. Modifications of the ND geometry and 
of dielectric properties of the material can lead to the real¬ 
ization of other domain configurations with topologically 
different hysteresis loops and even with a higher number 
of switching bits. 

To conclude, we performed the complete analysis of the 
domains formation and of the hysteresis switching pro¬ 
cess in the uniaxial ferroelectric ND as a function of the 
temperature and of the applied field. Various stable and 
metastable domain states were found. The most remark¬ 
able among them, the TBS, CCS and MDS, can be formed 
in a controlled way, during the regular zero-field cooling 
and systematic variation of the applied field. This paves 
the way to the predictable individual domain manipula¬ 
tion on the nanoscale. The discovered multibit switching 
opens new routes for the design of high-capacity nanosize 
memory-storage devices in the ferroelectric-based nano¬ 
electronics. 
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